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ABSTRACT 

Underbelly  blasts  to  vehicles  from  improvised  explosives  cause  severe  injuries  to  the  lower  extremities,  in¬ 
cluding  bone  fracture,  ligament  tear,  and  muscle  rupture.  While  these  injuries  may  seem  well-defined  through 
medical  imaging,  the  process  of  injury  and  the  effects  of  vehicle  system  design  to  protection  are  still  unclear.  In 
this  paper,  efforts  focused  on  developing  a  finite  element  model  of  the  lower  extremities  undergoing  high  strain 
rate  blast-induced  deformation  leading  to  injury  will  be  discussed.  A  hierarchical  approach  is  taken.  The 
process  of  high  strain  rate  direct  axial  loading  that  leads  to  bone  fracture  and  fragmentation  is  investigated. 

1.0  INTRODUCTION 

Improvised  explosive  devices  (IEDs)  pose  a  major  threat  to  the  Warfighter  and  are  oftentimes  hidden  along 
roads  and  detonated  underneath  trucks  or  tanks.  The  explosion  beneath  the  vehicle  is  commonly  referred  to 
as  an  underbelly  blast  event.  It  has  been  reported  that  8,159  IED  incidents  occurred  in  Afghanistan  in  2009, 
which  was  an  increase  from  3,867  and  2,677  which  were  reported  in  2008  and  2007,  respectively  [1].  During  an 
underbelly  blast,  the  lower  extremities  may  experience  accelerations  in  the  range  of  155  to  217  G  [2]  (note  that 
civilian  car  traffic  accidents  are  in  the  range  of  5  to  35  G  [3]).  In  order  to  design  enhanced  protective  systems 
for  this  particular  threat,  the  process  of  the  insult  to  the  injury  must  be  well  understood  to  discover  the  critical 
balance  between  performance  and  protection.  Computational  modeling  offers  a  powerful  tool  to  explore  the 
insult-to-injury  process  with  high-resolution. 

When  studying  a  complex  dynamic  process  such  as  this,  it  is  important  to  understand  the  length  and  time 
scales  with  the  event.  In  vehicular  mine  blasts,  multiple  types  of  blast  effects  are  present  including  primary, 
secondary  and  tertiary  blast  [4].  However,  many  of  the  musculoskeletal  injuries  are  currently  believed  to  be 
caused  by  tertiary  blast  effects  (for  example,  from  displacement  of  ah'  by  the  explosion  that  creates  a  blast 
wind  that  can  throw  victims  against  solid  objects)  [5].  During  an  underbelly  blast  event,  within  0.5  ms  of  the 
initiation  of  the  explosion,  a  shock  wave  hits  the  bottom  plate  of  the  vehicle  and  causes  a  rapid  rise  in  pressure. 
This  pressure  results  in  local  acceleration  and  subsequent  deformations  of  the  plate  which  apply  significant 
axial  loads  to  the  soldier  and  their  legs  [6].  As  the  blast  wave  reflects  under  the  vehicle,  a  pressure  force  acts  on 
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the  bottom  of  the  vehicle  and  causes  it  to  lift  off  the  ground.  The  acceleration  of  the  vehicle  and  the  collisions 
that  follow  are  also  a  significant  source  of  injury,  especially  if  the  soldier  is  not  properly  restrained.  The  body  is 
subjected  to  both  the  local  effects  (caused  by  shock  and  high  rate  deformation  of  the  vehicle)  and  global  effects 
(vehicle  motion  over  time)  as  a  result  of  the  mine  detonation  process.  The  nature  of  lower  extremity  injuries 
that  occur  during  a  mine  blast  may  be  broadly  distinguished,  similar  to  Kuppa  et  al.  [7],  as: 

•  Knee-thigh-hip  complex  fractures 

•  Knee  ligament  tears 

•  Tibia  plateau/condyle  fractures 

•  Tibia/fibula  shaft  fractures 

•  Calcaneus,  ankle,  and  midfoot  fractures 

•  Malleolar  fractures  and  ligament  tears  and  ruptures. 

Some  of  the  understanding  for  injuries  experienced  within  the  underbelly  blast  environment  can  be  obtained 
from  the  civilian  vehicle  crash  environment.  For  example,  fractures  of  the  ankle  and  lower  foot  which  are 
predominant  within  the  military  environment,  have  already  been  studied  in  the  past.  Manning  et  al.  [8]  report 
that  calcaneal,  talar,  midfoot  and  various  other  ankle  fractures  have  been  attributed  to  the  mechanism  of  axial 
loading  through  the  plantar  surface  of  the  foot.  While  these  injuries  might  be  expected  from  axial  loading, 
Crandall  et  al.  [9]  demonstrated  that  even  malleolar  fractures  can  occur  from  pure  dynamic  axial  loading  of  the 
leg.  Begemen  and  Paravasthu  [10]  also  conducted  experimental  impact  tests  where  unembalmed  cadaveric  legs 
were  subjected  to  a  uniaxial  plantar  surface  load  along  the  axis  of  the  tibia.  Pilon1  and  calcaneal  fractures  were 
observed  with  the  average  failure  tibial  axial  force  of  7590  N.  Kuppa  [7]  offers  a  useful  review  of  calcaneus, 
talus,  ankle  and  midfoot  fracture  research,  most  of  which  was  conducted  on  postmortem  human  subjects.  While 
this  past  research  highlights  some  common  features  with  military  blast  injuries,  it  is  important  to  note  that  the 
rate  of  loading  to  the  lower  extremities  may  be  10  -  100  times  greater  than  civilian  vehicular  crashes.  Thus, 
there  is  a  strong  possibility  that  different  injury  mechanisms  may  be  activated  during  these  events  and  new  or 
enhanced  injury  criteria  may  be  required. 

There  have  been  attempts  to  computationally  model  underbelly  blast  effects  in  the  lower  extremities.  Ni- 
lakantan  and  Tabiei  [2]  attempted  to  assess  injury  by  measuring  tibia  compressive  loads  following  an  under¬ 
belly  blast  event  using  finite  element  analysis  of  a  HYBRID  III  human  surrogate,  but  the  model  fails  to  capture 
specific  injury  mechanisms  such  as  bone  fracture  and  ligament  and  joint  failure.  They  make  clear  that  more  ad¬ 
vanced  and  biofidelic  lower  extremity  finite  element  models  are  needed.  Many  lower  extremity  studies  employ 
a  combined  experimental  and  simulation-based  approach.  For  example,  Manseau  and  Keown  [11]  compared  a 
physical  lower  leg  model  and  a  computer  model  in  the  development  of  enhanced  injury  assessment  methods, 
with  a  focus  on  the  ankle  complex  and  its  connections  to  the  tibia.  Once  again,  the  models  that  were  developed 
were  not  of  high  biological  fidelity  and  could  only  offer  some  information. 

In  this  paper  we  develop  a  series  of  numerical  models  in  an  effort  to  understand  dynamic  failure  of  the  lower 
extremities  subjected  to  underbelly  blast  loading.  The  lower  legs  will  be  modeled  using  three-dimensional, 
Lagrangian  finite  element  analysis  [12].  One  long-term  objective  of  this  research  is  to  have  the  capability  to 
provide  uncertainty  quantification  and  robust  validation  for  the  lower  extremity  model.  For  this  reason  we  adopt 
the  method  of  Thacker  et  al.  [13]  who  used  a  hierarchical  methodology  for  creating  a  numerical  model  of  the 
spine.  Section  2.0  describes  the  computational  methodology  and  results  for  three  distinct  finite  element  models: 
a  component  level  tibia  model  (2.2),  a  lower  leg  assembly  model  (2.3)  and  a  full  lower  extremity  system  model 
(2.4).  Conclusions  and  future  work  will  then  be  presented  in  Section  3.0. 

'A  pilon  fracture  occurs  when  one  bone  is  driven  into  another  bone  with  extreme  force  that  may  lead  to  comminuted  fracture.  The 
break  occurs  across  the  entire  bone  and  into  the  ankle  joint.  It  results  from  a  high-energy  loading  injury  from  the  foot  up  into  the  bone. 
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2.0  COMPUTATIONAL  METHODOLOGY  AND  RESULTS 
2.1  Model  Hierarchy 

Thacker  et  al.  [13]  provide  a  useful  hierarchical  verification  and  validation  scheme  that  decomposes  a  system 
level  model  into  components,  subassemblies  and  assemblies.  Within  the  validation  and  verification  process,  the 
process  should  start  at  smaller  length  scales,  such  as  the  components,  and  progress  to  the  last  stage  of  system 
level  models.  This  approach  requires  simulations  and  experiments  to  be  conducted  at  each  level.  While  the 
current  research  effort  adopts  the  approach,  we  note  that  this  effort  is  a  work-in-progress  and  is  still  in  an  early 
stage  of  building  such  a  infrastructure  and  research  team  to  support  both  the  experimental  and  numerical  models 
at  the  multiple  length  scales  required.  Nevertheless,  we  begin  to  define  the  subdivisions  for  the  lower  extremities 
for  this  effort  and  to  also  help  pave  the  way  for  future  work.  The  lower  extremity  system  is  partitioned  into 
components  (ligaments,  muscles,  bone,  skin),  subassemblies  (joints),  assemblies  (motion  segments  consisting 
of  muscles,  joints,  bones)  and  system  level  (whole  leg)  models  as  depicted  in  Figure  1. 


Real  World  Physical  System 
(Top-Level  of  Interest) 


Assemblies 


Subassemblies 


Components 


Figure  1 :  Validation  hierarchy  for  the  lower  extremities  based  on  the  method  of  Thacker  et  al.  [13].  The  system  is  decomposed 
into  components,  subassemblies,  assemblies  and  system  level  models. 


2.2  Component  Level  Modeling  of  Bone  Fracture  and  Failure  Mechanics 

Within  the  hierarchical  modeling  strategy,  the  component  is  the  smallest  length  scale  considered.  Here  we 
consider  only  one  component  of  many,  the  tibia  bone.  The  interest  in  bone  from  a  biomechanical  perspective 
is  largely  motivated  by  the  fact  that  bones  fracture,  and  that  different  fracture  mechanisms  may  be  activated 
during  a  blast  event,  as  described  in  Section  1 .0.  An  improved  understanding  of  the  mechanisms  underlying 
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bone  fracture  and  the  ability  to  predict  when  fractures  will  occur  may  help  in  the  development  of  new  strategies 
for  fracture  prevention  and  treatment.  Bones  can  fail  in  a  variety  of  ways.  If  a  bone  is  loaded  monotonically  to 
a  stress  level  that  exceeds  the  failure  or  ultimate  stress  of  the  bone  tissue,  then  fracture  will  occur.  This  type  of 
failure  can  occur  with  trauma  or  impact.  Bone  can  also  fail  at  stress  levels  well  below  the  yield  stress  through 
the  growth  of  preexisting  cracks,  but  this  will  not  be  discussed  here. 

In  this  component  level  study  the  tibia  is  segmented  into  two  parts  representing  cortical  and  trabecular  bone, 
with  the  properties  listed  in  Table  1.  The  medullary  cavity  and  other  anatomy  is  not  represented  here  as  the 
trabecular  bone  component  fills  the  entire  cortical  shell.  Axial  displacement  was  applied  to  the  proximal  (top) 
surface  of  the  tibia  as  the  distal  (bottom)  surface  remained  fixed. 

There  are  various  approaches  to  modeling  bone  fracture  numerically  including  element  deletion,  the  cohe¬ 
sive  zone  approach,  extended  Finite  Element  Modeling  (XFEM)  and  peridynamics.  As  Song  et  al.  [14]  point 
out,  the  element  deletion  method  is  one  of  the  simplest  methods  for  fracture  simulation  within  the  framework  of 
the  conventional  FEM  since  there  is  no  need  to  explicitly  represent  strong  discontinuities  in  displacement  fields 
(since  fracture  is  modeled  by  a  set  of  deleted  elements).  The  element  deletion  method  is  used  for  this  research 
due  to  simplicity.  Once  “broken”  the  elements  have  zero  stress  i.e.  zero  material  resistance,  which  is  imple¬ 
mented  by  using  a  constitutive  equation  in  which  the  stress  tends  to  zero  for  sufficiently  large  strain;  an  example 
of  such  a  stress-strain  law  is  shown  in  Figure  2.  Similar  to  other  numerical  failure  modeling  techniques,  element 
deletion  introduces  additional  length  and  time  scales  into  the  simulations  -  some  linked  with  physically-based 
parameters  such  as  material  toughness,  and  some  numerical,  such  as  the  softening  response.  Therefore,  effects 
of  constitutive  law  parameters  and  strain  rate  on  the  resulting  fracture  morphology  were  investigated. 


Figure  2:  Schematic  showing  the  elastic-fracture  constitutive  model  used  to  capture  fracture.  First  there  is  an  linear  elastic 
response  until  a  critical  stress,  <jc,  is  reached  within  the  element,  then  linear  stress  softening  occurs  until  the  element  critical 
crack  opening  stretch,  ec,  is  reached. 


2.2.1  Effect  of  Critical  Crack  Opening  Stretch 

One  of  the  crucial  points  in  using  the  element  deletion  method  is  the  scaling  of  the  constitutive  equation.  Unless 
the  constitutive  equation  is  adjusted  to  reflect  element  size,  the  released  energy  due  to  deleting  an  element 
depends  on  the  element  size,  which  causes  mesh  dependency.  To  reduce  mesh  dependency,  the  softening  curve 
slope  should  be  scaled  so  that  fracture  energy  is  independent  of  the  element  size.  The  energy  dissipation  in  an 
element  with  this  stress-strain  law  is  then  equated  to  the  surface  energy  of  a  crack  passing  through  the  element 
parallel  to  the  sides  by  modifying  the  stress-strain  law.  This  energy  consistency  renders  solutions  relatively 
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mesh  size  independent.  To  obtain  energy  equivalence  in  two  dimensional  problems  for  the  stress-strain  law 
shown  in  Figure  2,  the  upper  strain  limit  ec  (referred  to  as  the  critical  crack  opening  stretch)  is  scaled  so  that: 


Gche  =  l-Ee0ecAe 


(1) 


where  Gc  is  the  fracture  energy,  he  is  the  characteristic  dimension  (the  length  of  a  side  for  a  element)  and 
Ae  is  the  area  of  the  element  with  unit  thickness  assumed.  Recall  the  relationship  given  by  elastic  fracture 
mechanics  between  the  mode-I  stress  intensity  factor,  Kj,  and  the  energy  release  rate  of  the  material,  Gc  given 
by  Kj  =  E'GC  where  E'  =  E  for  plain  stress  or  FJ  =  E/(  1  —  v2)  for  plain  strain.  Ural  et  al.  [15]  reports 
that  bone  has  a  fracture  toughness  of  7-10  MPayTn.  Therefore,  in  plane  stress  with  E  =  19  GPa  and  Kjc  =  7 
MPa^/m,  Gc  =  2578  J/m2.  The  characteristic  dimension,  or  length  of  a  side  of  an  element  is  0.001  m  for  the 
tibia  finite  element  mesh.  In  order  to  calculate  an  estimate  for  the  upper  strain  limit  the  tibia  mesh  is  assumed  to 
consist  of  equilateral  tetrahedrons  in  which  the  area  of  the  elements  is  uniform  and  is  given  by  Ae  =  \  (0.001)2 
m2  =  5.0  x  10-7  m2. 


2 Gche  _  2(2578  J/m2) (0.001  m) 

acAe  ~  (100  MPa)(5.0  x  10-?m2) 


In  order  to  understand  the  effects  of  the  critical  crack  opening  stretch  on  fracture  morphology,  the  critical 
crack  opening  stretch  was  varied  between  0.1  and  0.001  m.  Figure  3  shows  mid-tibia  compressive  axial  stress 
plotted  in  terms  of  compressive  axial  strain  applied  to  the  ends  of  the  tibia  for  all  three  values  of  critical  crack 
opening  stretch  examined.  Once  the  crack  opening  critical  stretch  is  attained  the  elements  are  deleted  from  the 
computation,  thus  the  mass  of  the  elements  are  also  removed  from  the  total  system.  In  other  words,  the  total 
mass  of  the  system  decreases  as  fracture  increases.  Using  this  information,  a  simplistic  damage  parameter,  D, 
can  be  defined  as: 

D(t)  =  l-(  MfETEjM)  (3) 

\  initial  ) 

where  Mcurrent(t)  is  the  mass  of  the  system  at  time  t,  and  MimUai  is  the  mass  of  the  system  at  t  =  0.  The 
damage  defined  here  ranges  between  0  and  1,  with  1  representing  the  fully  damaged  system.  Damage  versus 
time  is  also  plotted  in  Figure  3. 

After  initial  compressive  loading,  stress  softening  occurs  at  0.2%  strain,  roughly  between  150-300  qs.  Note 
at  this  time  the  damage  parameter  also  begins  to  grow  indicating  crack  initiation.  Stress  values  for  the  initial 
softening  peak  are  similar  for  all  three  critical  crack  opening  stretch  values,  implying  that  the  fracture  initiation 
is  similar.  However,  secondary  stress  softening  that  occurs  between  0.4  and  0.6%  compressive  strain  arc  as 
much  as  26%  different,  with  noticeably  significant  values  for  damage.  This  stress  increase  and  then  softening 
continues  during  the  failure  process.  The  maximum  peak  compressive  axial  stress  reached  is  approximately  190 
MPa,  180  MPa  and  160  MPa  for  critical  crack  opening  stretch  values  of  0.1,  0.01  and  0.001,  respectively.  These 
stresses  correspond  to  nominal  axial  forces  of  approximately  41.4  kPa,  39.2  kPa  and  34.8  kPa,  respectively.  The 
rate  of  damage  growth  is  faster  for  smaller  values  of  critical  stretch  which  prevent  the  tibia  structure  to  support 
load  for  more  applied  strain  or  longer  duration  leading  to  higher  peak  loads  obtained. 

Figure  4  shows  fracture  profiles  for  all  three  values  of  critical  crack  opening  stretch.  Time  snapshots  of  the 
failure  process  is  also  shown.  For  the  case  of  ec  =  0.001  shown  in  Figure  4(a)  fracture  initiates  approximately 
0.483  ms,  however  that  fracture  is  near  the  location  of  applied  boundary  conditions.  At  1.861  ms  Figure  4(a) 
shows  how  the  fracture  has  initiated  near  the  mid-tibia  where  the  cross  sectional  area  is  minimum.  There 
are  three  microcracks  in  the  process  of  propagation.  At  2.344  ms  more  microcracks  are  initiated  and  are 
propagating,  and  some  microcracks  have  coalesed  with  others.  At  3.335  ms  almost  complete  coalescence  has 
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Figure  3:  Compressive  axial  stress  (solid)  and  damage  (dashed)  plotted  in  terms  of  axial  compressive  strain  and  time,  respec¬ 
tively,  for  various  critical  crack  opening  stretch  values. 


occurred  indicated  by  the  fragments  of  bone  isolated  in  space.  At  4.596  ms  the  bone  is  completely  fractured 
and  is  not  capable  of  carrying  any  load  due  to  the  fragmentation. 

A  similar  process  occurs  for  larger  values  of  ec,  however  since  the  critical  stretch  is  larger  crack  growth  is 
limited  which  thus  decreases  the  number  of  fragments  and  also  increases  fragment  size. 

2.2.2  Effect  of  Strain  Rate 

One  long-term  goal  of  this  research  is  to  accurately  predict  the  fracture  morphologies  of  the  bones  within  the 
lower  extremities  as  a  function  of  various  extrinsic  factors  such  as  loading  rate,  especially  since  the  mechanical 
properties  of  bone  have  been  reported  to  be  strain  rate  dependent.  Hansen  et  al.  [16]  measured  the  properties  of 
human  femoral  cortical  bone  at  strain  rates  ranging  between  0.14  -  29.1  s-1  in  compression  and  0.08  -  17  s~  1 
in  tension.  They  found  that  Young’s  modulus  generally  increased,  across  this  strain  rate  range,  for  both  tension 
and  compression. 

Figure  5  shows  the  compressive  axial  stress  versus  time  for  some  of  the  strain  rates  investigated  including 
10,  25  and  50  s^1.  The  strengthening  observed  in  these  simulations  is  similar  to  that  observed  in  experiments  on 
other  brittle  materials  and  is  believed  to  be  primarily  a  result  of  inertial  effects.  Unlike  some  micro-mechanical 
models,  inertial  effects  are  not  supplied  to  the  constitutive  description  ( e.g .  through  crack  growth  dynamics 
laws)  in  these  simulations.  Instead,  dilatation  of  the  highly  damaged  specimens  is  suppressed  because  material 
fragments  (which  have  mass)  resist  motion.  At  high  strain  rates  rapid  nucleation  and  coalescence  of  many 
cracks  is  a  result  of  large  stresses  that  develop  in  the  specimen  due  to  the  high  rate-of-deformation. 

As  shown  in  Figure  6,  the  damage  as  defined  by  Equation  3  has  an  increased  rate  of  growth  as  strain  rate 
increases.  The  total  amount  of  damage  also  increases  with  increasing  strain  rate.  Also  shown  in  the  figure 
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Figure  4:  Compressive  failure  response  for  the  tibia  subjected  to  compressive  loading  at  a  nominal  strain  rate  of  10  s-1  using  a 
critical  crack  opening  stretch  of  (a)  0.001,  (b)  0.01  and  (c)  0.1. 


are  snapshots  of  the  fragmentation  for  the  lowest  and  highest  strain  rate.  It  is  apparent  from  the  damage-time 
curves  that  crack  initiation  occurs  much  sooner  for  higher  rates  of  loading,  however  cracks  initiate  in  very 
different  spatial  locations.  For  lower  rates,  fracture  morphology  appears  to  more  transverse  near  the  mid-tibia 
whereas  for  higher  strain  rates  fracture  seems  to  more  proximal  and  could  possibly  more  described  as  a  pilon 
fracture.  The  transition  for  transverse  fracture  to  comminution  should  be  examined  further  to  gain  improved 
understanding  of  the  effects  of  bone  fragment  size  as  a  function  of  strain  rate.  From  qualitative  observations,  it 
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appears  that  perhaps  the  general  trend  is  that  as  the  strain  rate  increases  the  fragments  become  smaller.  In  this 
model,  the  smallest  allowable  fragment  size  is  the  elements.  Note  that  the  decrease  of  fragment  size  with  e  is 
also  predicted  by  a  number  of  theoretical  and  computational  methods  for  compressive  fragmentation  [17]. 


Figure  5:  Compressive  axial  stress  versus  time  for  strain  rates  of  10,  25  and  50  s-1. 


2.3  Modeling  the  Lower  Leg  Assembly 

In  accordance  with  our  hierarchical  validation  and  verification  strategy  presented  in  Section  2.1,  we  further 
consider  one  assembly  within  the  full  lower  extremity  model  the  lower  leg  of  the  region  below  the  knee.  Similar 
models  have  been  used  to  simulate  axial  impulsive  loading  applied  at  the  plantar  surface  of  the  foot  in  order  to 
understand  the  dynamic  response  of  the  ankle  related  to  injury  [18].  Following  Bandak  et  al.  [18]  the  model 
in  this  study  involves  the  simulation  of  cadaveric  lower  extremities  under  impulsive  axial  impact  that  attempt 
to  mimic  the  experiments  of  Yoganandan  et  al.  [19,  20].  In  these  experiments,  human  cadaveric  lower  legs 
with  the  ankle  initially  in  the  neutral  position  were  impacted  by  a  24.5  kg  pendulum  at  various  initial  velocities 
ranging  from  2.23  m/s  to  7.59  m/s  [19,  20].  A  similar  model  for  the  lower  leg  assembly  was  developed  for  this 
effort.  Figure  7a  shows  the  finite  element  model  of  the  lower  leg  assembly.  Figure  7b  shows  a  schematic  of  the 
computational  set-up  corresponding  to  the  experimental  model. 

For  the  lower  leg  assembly  model,  anatomy  distal  (below)  the  knee  is  represented.  Bones,  ligaments, 
muscles  and  skin  are  explicitly  represented  and  fully-segmented,  i.e.  each  bone,  muscle,  ligament,  etc.  can 
be  assigned  its  own  material  properties.  There  are  some  geometry  simplifications  that  were  taken  to  help 
aid  the  model  building  process.  In  reality,  bone  consists  of  many  layers  and  include  cortical  and  trabecular 
components,  plus  medullary  cavities  in  long  bones.  The  model  used  in  this  study  treats  bone  as  an  isotropic 
solid,  homogenous  body.  In  reality  the  foot  bone  structure  consists  of  26  bones,  with  numerous  articulations  and 
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Figure  6:  Damage  versus  time  for  various  strain  rates. 


ligaments.  The  model  used  in  this  study  uses  a  shrink-wrapped  mid-  and  forefoot  (the  talus,  navicular,  cuboid, 
cuneiforms,  metatarsals,  and  phalanges  are  fused).  The  calcaneus  was  not  altered  so  the  foot  is  represented  as 
two  objects:  “forefoot”  and  calcaneus.  In  reality,  13  muscles  run  along  the  tibia  and  fibula  in  the  lower  leg.  In 
the  model  only  4  of  those  muscles  are  explicitly  represented  currently;  the  tibialis  anterior,  tibialis  posterior, 
flexor  hallucis  longus  and  peroneus  brevis.  Tendons  are  not  identified  or  separated  from  muscle  bodies.  For 
the  other  soft  tissues,  a  bulk  material  was  utilized  to  encompass  the  remaining  muscles,  ligaments,  etc.  For 
example,  for  the  ankle  joints  gaps  between  bones  were  filled  in  with  the  bulk  material.  Articular  cartilage, 
joint  capsules,  synovial  fluid,  bursa,  etc.,  are  not  explicitly  represented.  All  components  of  the  model  (bones, 
muscles,  ligaments,  skin)  share  nodes  at  their  interfaces  i.e.  are  tied  together.  Finally,  a  rectangular  plate  was 
developed  and  placed  directly  underneath  and  in  contact  with  the  skin  of  the  foot.  This  plate  represented  the 
floor  of  the  truck  where  the  blast  would  occur  beneath.  During  the  simulations,  all  forces  affecting  the  leg 
model  would  first  travel  through  this  plate. 

Two  versions  of  the  lower  leg  model  with  different  bone  mesh  densities  (coarse  and  fine)  were  developed 
in  this  study.  A  comparison  between  the  two  can  be  visualized  in  Figure  8.  For  the  coarse  mesh,  there  were  a 
total  of  543,386  elements  in  the  model  which  includes  44,180  elements  for  the  tibia,  29,540  elements  for  the 
fibula,  32,522  elements  for  the  foot  and  5,713  elements  for  the  calcaneus.  For  the  fine  mesh,  there  were  a  total 
of  5,302,920  elements  in  the  model  which  includes  543,351  elements  for  the  tibia,  285,199  elements  for  the 
fibula,  83,233  elements  for  the  foot  and  54,370  elements  for  the  calcaneus.  The  simulations  of  the  lower  leg 
discussed  in  this  paper  were  performed  on  the  coarse  bone  mesh  model.  The  fine  mesh  model  will  be  used  for 
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Figure  7:  (a)  Lower  leg  model  geometry,  and  (b)  applied  velocity. 


future  study  for  better  visualization  of  crack  growth  and  propagation.  Similar  to  Bandak  et  al.  [18],  an  initial 
velocity  was  applied  to  the  rubber-plate  impactor.  In  this  study,  vD  =  10  m/s.  A  ballast  mass  of  11.3  kg  at  the 
top  (proximal)  end  of  the  tibia  was  used  in  accordance  with  Bandak  et  al.  [18]  to  simulate  the  weight  of  the 
experimental  slider  used  by  Yoganandan  et  al.  [19,  20]. 

Table  1  gives  a  list  of  the  material  constitutive  models  and  parameters  used  in  this  study.  Skeletal  muscle 
tissue  has  been  shown  to  exhibit  non-linear  anisotropic  viscoelastic  behavior  across  many  strain  rate  loading 
regimes.  In  this  model,  the  tissue  is  treated  as  non-linear  viscoelastic  using  a  viscoelastic  Swanson  model 
where  the  bulk  response  is  elastic  and  the  deviatoric  response  is  viscoelastic  [21,  22].  The  non-linear  elastic 
properties  are  obtained  from  the  recent  work  of  Lu  et  al.  who  represented  skeletal  muscle  tissue  as  a  composite 
of  an  isotropic  ground  substance  and  muscle  fibers  [23].  The  isotropic  matrix  was  modeled  using  a  hyperelastic 
neo-hookean  constitutive  law.  We  neglect  the  anisotropic  muscle  fiber  contribution  at  this  time.  Van  Loocke  et 
al.  examined  the  viscoelastic  properties  of  passive  skeletal  muscle  in  compression  [24]. 
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Figure  8:  Mesh  of  the  ankle  complex  within  the  lower  leg  assembly  model. 


Figure  9  shows  that  some  degree  of  load  sharing  occurs  between  the  tibia  and  fibula  when  the  lower  leg  is 
loaded  axially.  This  is  consistent  with  real-life  behavior  of  the  lower  leg,  and  is  a  result  unique  to  simulations 
of  actual  human  anatomy  (the  Hybrid-Ill  surrogate  and  its  corresponding  finite  element  models  do  not  have  a 
fibula).  In  the  time  series  contours,  it  can  be  seen  that  the  forefoot  experiences  the  majority  of  load  resulting  in 
significant  material  failure  and  element  deletion.  This  is  most  likely  an  artifact  of  the  simplified  geometry  used 
where  the  entire  forefoot  was  rigid;  none  of  the  articulations  between  the  metatarsals,  tarsals,  and  phalanges 
were  represented  in  this  model.  Furthermore  in  Figure  9,  it  can  be  seen  that  the  calcaneus  pivots  around  what 
would  be  the  ankle  joint  (displacements  are  magnified  by  a  factor  of  2  in  the  figure).  This  suggests  that  the  area 
representing  the  ankle  joint  in  this  model  must  be  stiffened  to  prevent  this  rotation  of  the  calcaneus  such  that 
crushing  fracture  on  the  planter  surface  of  the  heel  can  be  observed,  as  that  is  a  common  injury  experienced  in 
vehicular  underbelly  blasts  that  should  be  captured. 
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Table  1 :  Compilation  of  various  constitutive  models  and  parameters  used  for  the  lower  leg  simulation.  E  is  the  Young’s  Modulus, 
v  is  Poisson’s  ratio,  ac  is  the  critical  failure  stress,  cc  is  the  critical  crack  opening  stretch,  g^  is  the  long-term  shear  modulus, 
g„  are  prony  series  shear  terms,  p  is  the  density  and  K  is  the  bulk  modulus. 


Anatomic  Component 

Material  Model 

Material  Properties 

References 

Cortical  bone 

Elastic  with 

E  =  19.0  GPa 

[15] 

maximum 

< 

II 

O 

Id 

principle 

ac  =  132  MPa 

stress-based 

tc  —0.1 

Rone 

fracture  model 

p  =  1810  kg/m3 

Trabecular 

Elastic  with 

E  =  300  MPa 

[25] 

bone 

maximum 

v  =  0.45 

principle 

<rc=  1.5  MPa 

stress-based 

tc  —0.1 

fracture  model 

p  =  600  kg/m'5 

Isotropic  Visco- 

K  =  167  MPa 

[23,  26] 

Hyperelastic 

Cio  =  0.94  MPa  1 

(Mooney- 

Coi  =  0.94  MPa 

Rivlin) 

goo  =0.123 

[24] 

Muscles/Ligaments 

gi  =  0.465,  n  =  0.6  sec 
g2  =  0.2,  T2  =  6.0  sec 

[24] 

[24] 

g3  =  0.057,  T3  =  30.0  sec 

[24] 

g4  =  0.066,  T4  =  60.0  sec 

[24] 

g5  =  0.089,  Ts  =  300.0  sec 
p=  1000  kg/m3 

[24] 

Rubber  Shoe  Sole 

Elastic 

E  =  1  GPa 

v  =  0.4 

2.4  Real  World  Physical  System 

In  order  to  establish  a  system  of  models  ranging  from  individual  components  to  the  full  real  world  system, 
simplifications  were  made.  For  the  full  lower  extremity  model,  simple  geometry  consisting  of  cylinders  for  the 
bones  and  surrounding  muscle,  hexahedral  shapes  for  the  foot  and  non-friction  contact  interfaces  for  the  joints 
were  used  to  obtain  a  qualitative  response  force-time  curve  for  underbelly  loading.  Due  to  the  lack  of  anatomic 
geometry,  and  robust  constitutive  models  and  parameters,  this  approach  is  useful  to  understand  the  macroscopic 
impact  kinematics,  but  of  little  use  for  understanding  specific  injury  mechanisms.  Similar  to  the  computational 
study  of  Nilakantan  and  Tabiei  [2]  a  plate  which  represents  the  vehicle  floor  plate,  reaches  a  vertical  velocity  of 
10.7  m/s  after  5  ms.  While  Nilakantan  and  Tabiei  used  a  computational  model  of  a  Hybrid-Ill  surrogate,  here 
we  use  more  realistic  material  properties  for  bone  and  soft  tissue. 

As  shown  in  Figure  10,  within  approximately  2  ms  of  loading  the  axial  compressive  force  in  the  mid-tibia 
is  attained,  followed  by  a  sharp  rate  of  change  in  the  force-time  response  -  a  so-called  force  “plateau”.  This 
force  “plateau”  is  due  to  the  fact  that  the  velocity  of  the  tibia  has  reached  a  constant  value  because  the  upper 
leg  (between  the  knee  and  the  pelvis)  has  begun  to  travel  at  near  constant  velocity.  The  applied  loading  to  the 
floor  plate  ends  at  5  ms  at  which  time  the  force  response  diminishes  to  zero.  Figure  10  shows  results  of  this 
simplified  model  compared  to  the  dynamic  tibia  response  for  Hybrid-Ill  human  surrogates  [2]  and  cadaveric 
experiments  [27]  using  similar  boundary  conditions. 

When  including  muscle  components  into  the  skeletal  model,  the  materials  shear  modulus  became  an  impor¬ 
tant  contributor  to  the  axial  compressive  force  inflicted  on  the  tibia.  These  experiments  showed  that  the  shear 
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Figure  9:  Simulation  contours  of  the  lower  leg  (course  model)  subjected  to  vertical  loading  process  of  10  s~L 


modulus  within  the  muscle  could  control  behavior  within  the  initial  1  ms  of  the  tibial  injury  event.  As  shown  in 
Figure  1 1  high  shear  modulus  for  muscle  was  made  an  upper  bound  by  making  it  26.8  MPa.  Low  shear  modulus 
was  controlled  by  a  lower  bound  of  26.8  kPa.  The  magnitude  of  compressive  force  within  this  time  frame  was 
maximal,  around  2000  N,  when  the  muscle  exhibited  high  shear  properties.  When  the  muscle  exhibited  low 
shear  muscle,  or  no  muscle  at  all,  the  force  was  half  as  much  as  its  high  shear  counterpart.  It  was  found  that 
this  correlation  was  indirectly  related  to  the  muscles  effect  on  rate  of  accelerative  loading.  High-shear  muscle 
led  to  a  slower  rate  of  tibial  acceleration,  due  to  the  muscle’s  inertial  effects.  The  low-shear  and  no  muscle’s 
rate  of  acceleration  was  nearly  three  times  that  of  high-shear  muscle.  This  result  shows  that  there  is  a  clear  link 
between  rate  of  accelerative  loading  and  compressive  force  levels  within  the  tibia. 

While  this  approach  fails  to  give  insight  into  the  mechanics  of  any  specific  lower  extremity  injury,  it  does 
provide  insight  into  the  bio-kinematics  of  the  lower  extremities  during  an  explosive  event.  The  approach  was 
also  very  useful  to  understand  the  importance  of  the  contributions  of  individual  components  such  as  the  muscle 
or  joints  to  the  overall  real  world  physical  system  model.  In  addition,  a  simplified  model  that  offers  macroscopic 
results  quickly,  such  as  this,  could  be  use  to  develop  boundary  conditions  for  the  smaller  length  scale  models 
including  assemblies  and  subassemblies.  For  example,  it  is  important  to  note  that  the  mid-tibia  axial  force 
reaches  a  maximum  value  within  2  ms  of  loading. 
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Figure  10:  For  the  full  lower  extremity  model,  simple  geometry  consisting  of  cylinders  for  the  bones  and  surrounding  muscle, 
hexahedral  shapes  for  the  foot  and  non-friction  contact  interfaces  for  the  joints  were  used  to  obtain  a  qualitative  response  force¬ 
time  curve  for  underbelly  loading.  Results  are  compared  to  simulations  of  Hybrid-Ill  surrogates  [2]  and  cadaveric  experiments 
[27]. 


3.0  CONCLUSIONS  AND  FUTURE  WORK 

This  paper  has  described  a  research  effort  to  develop  a  hierarchical  modeling  approach  for  the  lower  extremities 
subjected  to  military  underbelly  blasts  and  has  specifically  focused  on  developing  a  finite  element  model  of  the 
lower  extremities  undergoing  high  strain  rate  blast-induced  deformation.  At  the  component  level,  the  process 
of  bone  fracture  and  fragmentation  from  high  strain  rate  axial  loading  to  the  tibia  was  investigated.  The  effect 
of  critical  crack  opening  stretch  and  strain  on  the  resulting  leg  fracture  was  studied  and  had  significant  influence 
on  fracture  morphology.  Apparent  strengthening  was  observed  for  higher  rates  of  loading,  which  are  attributed 
to  inertial  effects.  At  the  assembly  level,  a  lower  leg  consisting  of  many  different  tissues  was  developed  and 
evaluated.  Similar  to  past  experimental  research,  we  numerically  observe  calcaneal,  talar,  midfoot  fractures, 
as  well  as  ankle  rotation  before  any  tibia  or  fibula  fracture  is  observed.  At  the  system  (full  lower  extremity) 
level,  simulations  consisting  of  simplified  geometry  with  more  realistic  material  constitutive  laws  as  compared 
to  Hybrid  III  finite  element  models  were  developed.  The  models  captured  the  kinematics  of  the  lower  leg  and 
the  axial  tibia  force  was  compared  to  cadaveric  experimental  tests  as  well  as  results  from  a  Hybrid  III  human 
surrogate  finite  element  model.  All  three  models  need  to  be  studied  in  more  depth  and  refined.  For  example, 
for  each  component  within  the  lower  leg  assembly  a  computational  model  should  be  developed  and  analyzed  - 
both  computationally  and  experimentally.  To  summarize,  the  following  offers  a  list  of  future  work  that  should 
be  completed: 

•  Experimental  and  computational  component  level  studies  including  bone,  ligaments,  muscle  and  skin  at 
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0.002 


0.004  0.006 

Time  (ms) 


0.008 


0.01 


Figure  11 :  Effect  of  soft  tissue  shear  modulus  on  the  axial  compressive  force. 


the  appropriate  strain  rates. 

•  Include  more  complex  material  models  using  anisotropic  and  poro- visco-elasticity  that  capture  the  real¬ 
istic  tissue  response. 

•  Examine  additional  assemblies,  such  as  the  knee-thigh-pelvis  complex. 

•  Develop  computational  model  of  lower  extremities  that  is  linked  with  the  spine  to  elucidate  injury  mech¬ 
anism  associated  with  coupled  loading. 

•  Improve  the  computational  scalability  and  performance  on  supercomputing  resources. 

•  Improve  availability  of  human  data  pertaining  to  lower  leg  impact  during  a  mine  blast  underneath  armored 
vehicles,  especially  testing  on  human  cadavers  in  order  to  better  understand  the  injury  assessment  and 
establish  a  reliable  and  extensive  source  of  data. 

•  Leverage  detailed  understanding  of  injury  criteria  obtained  from  computational  studies  to  further  refine 
military  blast  injury  criteria. 

•  Use  injury  mechanism-based  understanding  to  improve  protection  system  design. 

In  conclusion,  computational  failure  modeling  of  the  lower  extremities  offers  valuable  opportunities  to  improve 
protection  for  the  Warfighter. 
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